home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / min / test.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-11-01  |  6.3 KB  |  216 lines
  1. /* min/test.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19. #include <config.h>
  20. #include <stdlib.h>
  21. #include <gsl/gsl_math.h>
  22. #include <gsl/gsl_min.h>
  23. #include <gsl/gsl_errno.h>
  24. #include <gsl/gsl_test.h>
  25. #include <gsl/gsl_ieee_utils.h>
  26.  
  27. #include "test.h"
  28. #include "min.h"
  29.  
  30. /* stopping parameters */
  31.  
  32. const double EPSABS = 0.001 ;
  33. const double EPSREL = 0.001 ;
  34.  
  35. const unsigned int MAX_ITERATIONS = 100;
  36.  
  37. void my_error_handler (const char *reason, const char *file,
  38.                int line, int err);
  39.  
  40. #define WITHIN_TOL(a, b, epsrel, epsabs) \
  41.  (fabs((a) - (b)) < (epsrel) * GSL_MIN(fabs(a), fabs(b)) + (epsabs))
  42.  
  43. int
  44. main (void)
  45. {
  46.   gsl_function F_cos, F_func1, F_func2, F_func3, F_func4;
  47.   
  48.   const gsl_min_fminimizer_type * fminimizer[4] ;
  49.   const gsl_min_fminimizer_type ** T;
  50.  
  51.   gsl_ieee_env_setup ();
  52.  
  53.   fminimizer[0] = gsl_min_fminimizer_goldensection;
  54.   fminimizer[1] = gsl_min_fminimizer_brent;
  55.   fminimizer[2] = 0;
  56.  
  57.   F_cos = create_function (f_cos) ;
  58.   F_func1 = create_function (func1) ;
  59.   F_func2 = create_function (func2) ;
  60.   F_func3 = create_function (func3) ;
  61.   F_func4 = create_function (func4) ;
  62.  
  63.   gsl_set_error_handler (&my_error_handler);
  64.  
  65.   for (T = fminimizer ; *T != 0 ; T++)
  66.     {
  67.       test_f (*T, "cos(x) [0 (3) 6]", &F_cos, 0.0, 3.0, 6.0, M_PI);
  68.       test_f (*T, "x^4 - 1 [-3 (-1) 17]", &F_func1, -3.0, -1.0, 17.0, 0.0);
  69.       test_f (*T, "sqrt(|x|) [-2 (-1) 1.5]", &F_func2, -2.0, -1.0, 1.5, 0.0);
  70.       test_f (*T, "func3(x) [-2 (3) 4]", &F_func3, -2.0, 3.0, 4.0, 1.0);
  71.       test_f (*T, "func4(x) [0 (0.782) 1]", &F_func4, 0, 0.782, 1.0, 0.8);
  72.  
  73.       test_f_e (*T, "invalid range check [4, 0]", &F_cos, 4.0, 3.0, 0.0, M_PI);
  74.       test_f_e (*T, "invalid range check [1, 1]", &F_cos, 1.0, 1.0, 1.0, M_PI);
  75.       test_f_e (*T, "invalid range check [-1, 1]", &F_cos, -1.0, 0.0, 1.0, M_PI);
  76.     }
  77.   test_bracket("cos(x) [1,2]",&F_cos,1.0,2.0,15);
  78.   test_bracket("sqrt(|x|) [-1,0]",&F_func2,-1.0,0.0,15);
  79.   test_bracket("sqrt(|x|) [-1,-0.6]",&F_func2,-1.0,-0.6,15);
  80.   test_bracket("sqrt(|x|) [-1,1]",&F_func2,-1.0,1.0,15);
  81.  
  82.   exit (gsl_test_summary ());
  83. }
  84.  
  85. void
  86. test_f (const gsl_min_fminimizer_type * T, 
  87.         const char * description, gsl_function *f,
  88.     double lower_bound, double middle, double upper_bound, 
  89.         double correct_minimum)
  90. {
  91.   int status;
  92.   size_t iterations = 0;
  93.   double m, a, b;
  94.   double x_lower, x_upper;
  95.   gsl_min_fminimizer * s;
  96.  
  97.   x_lower = lower_bound;
  98.   x_upper = upper_bound;
  99.  
  100.   s = gsl_min_fminimizer_alloc (T) ;
  101.   gsl_min_fminimizer_set (s, f, middle, x_lower, x_upper) ;
  102.   
  103.   do 
  104.     {
  105.       iterations++ ;
  106.  
  107.       status = gsl_min_fminimizer_iterate (s);
  108.  
  109.       m = gsl_min_fminimizer_minimum(s);
  110.       a = gsl_min_fminimizer_x_lower(s);
  111.       b = gsl_min_fminimizer_x_upper(s);
  112.  
  113. #ifdef DEBUG
  114.       printf("%.12f %.18f %.12f %.18f %.12f %.18f\n", 
  115.              a, GSL_FN_EVAL(f, a), m, GSL_FN_EVAL(f, m), b, GSL_FN_EVAL(f, b));
  116. #endif
  117.  
  118.       if (a > b)
  119.     gsl_test (GSL_FAILURE, "interval is invalid (%g,%g)", a, b);
  120.  
  121.       if (m < a || m > b)
  122.     gsl_test (GSL_FAILURE, "m lies outside interval %g (%g,%g)", m, a, b);
  123.  
  124.       if (status) break ;
  125.  
  126.       status = gsl_min_test_interval (a, b, EPSABS, EPSREL);
  127.     }
  128.   while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS);
  129.  
  130.   gsl_test (status, "%s, %s (%g obs vs %g expected) ", 
  131.         gsl_min_fminimizer_name(s), description, 
  132.         gsl_min_fminimizer_minimum(s), correct_minimum);
  133.  
  134.   /* check the validity of the returned result */
  135.  
  136.   if (!WITHIN_TOL (m, correct_minimum, EPSREL, EPSABS))
  137.     {
  138.       gsl_test (GSL_FAILURE, "incorrect precision (%g obs vs %g expected)", 
  139.         m, correct_minimum);
  140.     }
  141.  
  142.   gsl_min_fminimizer_free (s);
  143.  
  144. }
  145.  
  146. void
  147. test_f_e (const gsl_min_fminimizer_type * T, 
  148.       const char * description, gsl_function *f,
  149.       double lower_bound, double middle, double upper_bound, 
  150.           double correct_minimum)
  151. {
  152.   int status;
  153.   size_t iterations = 0;
  154.   double x_lower, x_upper;
  155.   double a, b;
  156.   gsl_min_fminimizer * s;
  157.  
  158.   x_lower = lower_bound;
  159.   x_upper = upper_bound;
  160.  
  161.   s = gsl_min_fminimizer_alloc (T) ;
  162.   status = gsl_min_fminimizer_set (s, f, middle, x_lower, x_upper) ; 
  163.  
  164.   if (status != GSL_SUCCESS) 
  165.     {
  166.       gsl_min_fminimizer_free (s) ;
  167.       gsl_test (status == GSL_SUCCESS, "%s, %s", T->name, description);
  168.       return ;
  169.     }
  170.  
  171.   do 
  172.     {
  173.       iterations++ ;
  174.       gsl_min_fminimizer_iterate (s);
  175.       a = gsl_min_fminimizer_x_lower(s);
  176.       b = gsl_min_fminimizer_x_upper(s);
  177.  
  178.       status = gsl_min_test_interval (a, b, EPSABS, EPSREL);
  179.     }
  180.   while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS);
  181.  
  182.   gsl_test (!status, "%s, %s", gsl_min_fminimizer_name(s), description, 
  183.         gsl_min_fminimizer_minimum(s) - correct_minimum);
  184.  
  185.   gsl_min_fminimizer_free (s);
  186. }
  187.  
  188. void
  189. my_error_handler (const char *reason, const char *file, int line, int err)
  190. {
  191.   if (0)
  192.     printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err);
  193. }
  194.  
  195. int
  196. test_bracket (const char * description,gsl_function *f,double lower_bound, 
  197.           double upper_bound, unsigned int max)
  198. {
  199.   int status;
  200.   double x_lower, x_upper;
  201.   double f_upper,f_lower,f_minimum;
  202.   double minimum;
  203.  
  204.   x_lower=lower_bound;
  205.   x_upper=upper_bound;
  206.   SAFE_FUNC_CALL (f,x_lower,&f_lower);
  207.   SAFE_FUNC_CALL (f,x_upper,&f_upper);
  208.   status=gsl_min_find_bracket(f,&minimum,&f_minimum,&x_lower,&f_lower,&x_upper,&f_upper,max);
  209.   gsl_test (status,"%s, interval: [%g,%g], values: (%g,%g), minimum at: %g, value: %g",
  210.         description,x_lower,x_upper,f_lower,f_upper,minimum,f_minimum);
  211.   return status;
  212. }
  213.  
  214.  
  215.  
  216.